home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Pascal Super Library
/
Pascal Super Library (CW International)(1997).bin
/
MATH
/
PMAT12
/
PMAT.EXE
/
TESTOP.PAS
< prev
Wrap
Pascal/Delphi Source File
|
1993-01-17
|
8KB
|
280 lines
Program testop;
{ test suite for pmatop }
Uses pmat, pmatop;
Procedure testMfunctions;
Var
x,y : vmatrixptr;
Begin
Dispatch^.Inclevel;
new( x, makematrix( 1, 1 ) );
new( y, makematrix( 1, 1 ) );
x := matequals( x, add( Ident( 5 ), fill( 5, 5, 1 ) ) );
{ sort rows in ascending order of the first column }
y := matequals( y, MSort( x, 1, 1 ) );
y^.Show( 'sorted on rows' );
{ sort columns in ascending order of the second row }
y := matequals( y, MSort( y, 2, 0 ) );
y^.Show( 'sorted on columns' );
y := matequals( y, Mexp( x ) );
y^.Show( 'exponent of elements' );
y := matequals( y, Mabs( x ) );
y^.Show( 'absolute value of elements' );
y := matequals( y, Mlog( x ) );
y^.Show( 'natural logrithm of elements' );
y := matequals( y, Msqrt( x ) );
y^.Show( 'square root of elements' );
dispose( x, killvmatrix );
dispose( y, killvmatrix );
Dispatch^.Declevel;
End;
Procedure testPatterned;
Var
H, K : vmatrixptr;
Begin
Dispatch^.Inclevel;
new( H, makematrix( 5, 5 ) );
new( K, makematrix( 5, 5 ) );
h := matequals( h, helm( 5 ) );
h^.Show( 'Helmert Matrix' );
writeln;
writeln( 'trace of the helmert(5): ', trace( H ): 7: 5 );
writeln( 'Determinant of Helmert Matrix: ', Det( H ): 7: 5 );
k := matequals( k, index( 1, 5, 1 ) );
k^.Show( '1 5 1' );
k := matequals( k, index( 5, 1, - 1 ) );
k^.Show( '5 1 -1' );
k := matequals( k, index( 1, 5, 2 ) );
k^.Show( '1 5 2' );
k := matequals( k, index( 5, 1, - 2 ) );
k^.Show( '5 1 -2' );
k := matequals( k, Kron( Ident( 2 ), H ) );
k^.Show( 'I@H' );
k := matequals( k, H );
k^.mm( 1, 1 )^ := 1.0E - 9;
writeln;
writeln( 'Determinant of Pivoted Matrix: ', Det( k ): 7: 5 );
dispose( h, killvmatrix );
dispose( k, killvmatrix );
Dispatch^.Declevel;
End;
Procedure testDecomp;
Var
H, S, Q, D, R : vmatrixptr;
Begin
Dispatch^.Inclevel;
new( H, makematrix( 5, 5 ) );
new( S, makematrix( 5, 5 ) );
new( Q, makematrix( 5, 5 ) );
new( D, makematrix( 5, 5 ) );
new( R, makematrix( 5, 5 ) );
h := matequals( h, Add( Ident( 5 ), Fill( 5, 5, 1 ) ) );
S := matequals( S, Cholesky( h ) );
S^.Show( 'Cholesky decomposition' );
S := matequals( S, mult( Tran( S ), S ) );
S^.Show( 'S`S' );
QR( h, Q, R, true );
Q^.Show( 'Q' );
R^.Show( 'R' );
S := matequals( S, mult( Q, R ) );
S^.Show( 'QR' );
SVD( H, Q, D, R, TRUE, TRUE );
Q^.Show( 'S' );
D^.Show( 'D' );
R^.Show( 'R' );
S := matequals( S, mult( mult( Q, Diag( D ) ), Tran( R ) ) );
S^.Show( 'SDR`' );
S := matequals( S, Ginv( H ) );
S^.Show( 'Ginv(H)' );
S := matequals( S, mult( S, H ) );
S^.Show( 'h*Ginv(h)' );
dispose( h, killvmatrix );
dispose( S, killvmatrix );
dispose( Q, killvmatrix );
dispose( D, killvmatrix );
dispose( R, killvmatrix );
Dispatch^.Declevel;
End;
Procedure testShape;
Var
H, S, Q, D, R : vmatrixptr;
Begin
Dispatch^.Inclevel;
new( H, makematrix( 1, 1 ) );
new( S, makematrix( 1, 1 ) );
h := matequals( h, Add( Ident( 5 ), Fill( 5, 5, 1 ) ) );
s := matequals( s, Vec( h ) );
S^.Show( 'Vec(h)' );
s := matequals( s, shape( s, 5 ) );
S^.Show( 'Shape(Vec(s)),5' );
s := matequals( s, Diag( h ) );
s^.show( 'Diagonal elements of h' );
Dispatch^.Declevel;
End;
Procedure testFFT;
Var
i,n : integer;
omega, sigma, xbar : double;
test,f,t, average, centered, power_spectrum : vmatrixptr;
covar, periodogram : vmatrixptr;
Begin
n := 30; { n=2*3*5 }
omega := pi / n;
new( t, makematrix( n, 1 ) );
new( f, makematrix( n, 2 ) );
new( test, makematrix( n, 2 ) );
new( average, makematrix( 1, 1 ) );
new( centered, makematrix( 1, 1 ) );
new( power_spectrum, makematrix( 2 * n, 1 ) );
new( covar, makematrix( 2 * n, 1 ) );
new( periodogram, makematrix( n, 1 ) );
For i := 0 To n - 1 Do
t^.mm( i + 1, 1 )^ := sin( i * omega ) + cos( 0.3 * i * omega );
f := matequals( f, fft( t, TRUE ) );
f^.show( 'time to frequency transform' );
test := matequals( test, ch( t, fft( f, FALSE ) ) );
test^.show( 'frequency to time transform' );
{ get power spectrum, and serial correlations }
average := matequals( average, mult( Fill( 1, t^.r, 1 ), t ) );
xbar := average^.m( 1, 1 ) / n;
centered := matequals( centered, sub( t, Fill( t^.r, 1, xbar ) ) );
centered := matequals( centered, cv( t, Fill( t^.r, 1, 0 ) ) );
f := matequals( f, fft( centered, TRUE ) );
For i := 1 To 2 * n Do
power_spectrum^.mm( i, 1 )^ := (sqr( f^.m( i, 1 ) ) + sqr( f^.m( i, 2 ) ) ) / n;
power_spectrum^.show( 'power spectrum' );
covar := matequals( covar, fft( power_spectrum, FALSE ) );
sigma := covar^.m( 1, 1 );
covar := matequals( covar, submat( covar, 1, n, 1, 1 ) );
For i := 1 To n Do Begin
covar^.mm( i, 1 )^ := covar^.m( i, 1 ) / sigma;
periodogram^.mm( i, 1 )^ := power_spectrum^.m( i, 1 ) / sigma;
End;
covar^.show( 'serial correlations' );
periodogram := matequals( periodogram, Mlog( periodogram ) );
periodogram^.show( 'periodogram' );
dispose( t, killvmatrix );
dispose( f, killvmatrix );
dispose( test, killvmatrix );
dispose( average, killvmatrix );
dispose( centered, killvmatrix );
dispose( power_spectrum, killvmatrix );
dispose( covar, killvmatrix );
dispose( periodogram, killvmatrix );
End;
Procedure TestSums;
Var
a, b : vmatrixptr;
Begin
new( a, makematrix( 1, 1 ) );
new( b, makematrix( 1, 1 ) );
a := matequals( a, Fill( 5, 5, 1 ) );
b := matequals( b, Sum( a, ALL ) );
b^.show( 'sum over all' );
b := matequals( b, Sum( a, ROWS ) );
b^.show( 'sum over rows' );
b := matequals( b, Sum( a, COLUMNS ) );
b^.show( 'sum over columns' );
b := matequals( b, Sumsq( a, ALL ) );
b^.show( 'sum sq over all' );
b := matequals( b, Sumsq( a, ROWS ) );
b^.show( 'sum sq over rows' );
b := matequals( b, Sumsq( a, COLUMNS ) );
b^.show( 'sum sq over columns' );
b := matequals( b, Cusum( a ) );
b^.show( 'Cusum of a' );
a := matequals( a, b );
b := matequals( b, Mmax( a, ALL ) );
b^.show( 'max over all' );
b := matequals( b, Mmax( a, ROWS ) );
b^.show( 'max over rows' );
b := matequals( b, Mmax( a, COLUMNS ) );
b^.show( 'max over columns' );
b := matequals( b, Mmin( a, ALL ) );
b^.show( 'min over all' );
b := matequals( b, Mmin( a, ROWS ) );
b^.show( 'min over rows' );
b := matequals( b, Mmin( a, COLUMNS ) );
b^.show( 'min over columns' );
dispose( a, killvmatrix );
dispose( b, killvmatrix );
End;
Procedure testElementary;
Var
a : vmatrixptr;
Begin
new( a, makematrix( 1, 1 ) );
a := matequals( a, Add( Ident( 5 ), Fill( 5, 5, 1 ) ) );
Crow( a, 2, 3.0 );
a^.show( 'r2*3.0' );
Srow( a, 1, 2 );
a^.show( 'r1 <-> r2' );
Lrow( a, 1, 2, - 3.0 );
a^.show( '-3.0*r2 + r1' );
Ccol( a, 2, 3.0 );
a^.show( 'c2*3.0' );
Scol( a, 1, 2 );
a^.show( 'c1 <-> c2' );
Lcol( a, 1, 2, - 3.0 );
a^.show( '-3.0*c2 + c1' );
dispose( a, killvmatrix );
End;
Begin
{main function}
testMfunctions;
testPatterned;
testDecomp;
testShape;
testFFT;
testSums;
testElementary;
End.